Submitted to Annals of Applied Probability 



GLOBAL REGULARITY AND PROBABILISTIC SCHEMES 
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rvj ' In a rather general setting of multivariate stochastic volatility 

market models we derive global iterative probabilistic schemes for 
computing the free boundary and its Greeks for a generic class of 
American derivative models using front-fixing methods. Establish- 

J~^ , ment of convergence is closely linked to a proof of global regularity 

of the free boundary surface. 

< 

1^ ' 1. Introduction. The need of probabilistic representations of solutions 

,^ ■ of problems in finance is closely linked to the curse of dimensionality which 

C^ , enforces the application of Monte-Carlo methods especially for more intri- 

cate problems. The situation is especially challenging for nonlinear problems 
such as the free boundary problems related to pricing, hedging, and the op- 
timal exercise strategies of American options. The present work establishes 
^ I a stable iterative scheme which generates a sequence functions and their 

derivatives up to second order converging to the free boundary function 

f^ . (the image of which is the free boundary surface) and its derivatives for a 

standard class of American derivatives modelled in a quite general frame- 

t^ ' work of diffusion type market models. Hence, it contributes to all these 

three financial problems. However, it is also a contribution to the theoret- 

f^ i ical question of regularity of the free boundary itself where we establish 

results which are neither covered by the general results of [">, (>] (cf. dis- 
cusson below) nor by more special recent regularity results on American 
derivatives in [25, 28]. There are several methods in order to investigate 
C^ I regularity used in the literature. One is to differentiate equations in order 

to get equations for derivatives. Another is to find estimates for tangent 
paraboloids of the solutions (cf. [30, 31, 32]. Another is to find the solution 
in explicit form which is possible only in a very limited class of problems. 
In [6] the multivariate case of parabolic potential problem is investigated, 
i.e. the case of constant coefficients and with no assumption on the sign of 
the solution. The points of the free boundary are classified into regular and 
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2 J. KAMPEN 

singular points by an energy and a density criterion and C°° regularity is 
proved around the regular points. Singular points are not considered. Hence, 
the regularity results of the present work which establishes overall regularity 
especially for the American basket Put option are not covered even in the 
case of the multivariate Black-Scholes model. However our results extend to 
a considerable class of models with variable coefficients. The only condition 
(beside regularity conditions on the coefficients) is a global graph condition 
which allows for a global transformation. Our method is based on the front 
fixing method. The term front-fixing refers to a nonlinear transformation 
to fix the free boundary and to solve the resulting non-linear problem. The 
transformation depends essentially on the data (in financial models on the 
final data) . The method of front fixing itself is discussed in several papers on 
univariate American options (cf. [23, 33]. In the univariate case these meth- 
ods are less efficient than a moving boundary approach (cf. [22]). We shall 
investigate the American basket put option as the most prominent example 
treated in the literature. Regularity results for the American put option in 
the univariate case go back to [20] and [29] in the case of constant coef- 
ficients. In [11] the univariate case where coefficients depend on time and 
space. For a more recent contribution in the univariate case with variable 
coefficients cf. [3]. In the present work we obtain global regularity results 
for the American basket put. However, our method can be applied to the 
other most popular options (we shall provide a list below). In case of the 
Basket Put option we transform the free boundary problem to a nonlinear 
problem on a domain homeomorphic to a cube. We present the solution of 
the free boundary surface in terms of a nonlinear integral equations involv- 
ing convolutions with transition densities of linear parabolic equations and 
linear Volterra equations. Analysis of this equations leads to regularity re- 
sults which are simple (use of Banach fixed point theorem) and charming 
except that they involve some use of the Levy expansion of linear parabolic 
equations. Our main result here is the global regularity of the free boundary 
surface except for the final time T if a global graph condition is satisfied 
(which can be justified for a considerable class of models). The nonlinear in- 
tegral equation is also the basis for our probabilistic scheme for computing 
the free boundary function, its time derivative and its spatial derivatives up 
to second order (Greeks). The only other Monte-Carlo method for comput- 
ing the free boundary (not Greeks) of a multivariate American Put option 
known to the author is [2, 21]. We break down the solution of this inte- 
gral equation to an iteration of the solution of linear equations (essentially 
two parabolic equations and a linear Volterra integral equation). The cor- 
responding convolutions with the transition density can be computed by 
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FREE BOUNDARIES AND THEIR GREEKS 3 

WKB approximations (cf. [16, 17, 18]). Recent developments in the compu- 
tations of Greeks with Monte-Carlo Methods (cf. [7, 13, 15]. Especially in 
[15] it is shown that WKB approximations can be a very efficient tool for 
the computation of Greeks for high-dimensional models. The convergence of 
the resulting iteration method is based on the regularity results for the free 
boundary. 

The outline of the present article is as follows. In Section 2 we describe 
the general frame work of market models considered. We restrict ourselves 
to diffusion models (possibly with stochastic volatility). However extension 
especially to Levy models is possible. We formulate the free boundary prob- 
lem to be solved. However the treatment of global operators in the context 
of free boundaries leads to additional involvements which will be consid- 
ered elsewhere. In Section 3 we apply the front fixing method to the market 
models set up in the previous Section. We apply it especially for the Amer- 
ican Put option (the most popular payoff function in the literature), but 
mention also other standard payoffs to which it may be applied. The free 
boundary problem set up in the previous section is globally transformed to 
a nonlinear problem on a hypercube. The free boundary function appears 
in the coefficients of this nonlinear parabolic function and is coupled to 
the boundary conditions. Especially the smooth fit condition of the origi- 
nal free boundary transforms is code into a mixed boundary condition on a 
hyperplane. The involved global graph condition which is formulated. Note 
that the possible failure of the smooth fit condition for general Levy mod- 
els and beyond is one of the mentioned difficulties beyond market models 
with path-continuous models which make a detailed separated treatment of 
these models necessary. In Section 4 we derive a nonlinear integral equation 
which characterizes the free boundary and is based on the nonlinear initial 
boundary value problem characterizing the value function and established 
in the previous Section. In Section 5 a proof of the global existence and the 
regularity of the free boundary function is established. This is done via the 
nonlinear integral equation established in the previous Section. The integral 
equation is shown to determine a map of the Banach space of Holder contin- 
uous functions with Holder continuous derivatives up to the second order. 
Note that interior norms are used because the problem is not differentiable 
at maturity time and only C^ on the boundary hyperplane where the mixed 
boundary condition is related to smooth fit condition of the original prob- 
lem. The solution of the free boundary function is the fixed point of the 
nonlinear integral equation in the indicated Banach space. Global existence 
follows easily by standard theory of linear parabolic problems. In Section 6 
an iteration scheme for the computation of the free boundary function and 

imsart-aap ver. 2007/12/10 file: aapamerican.tex date: July 7, 2008 



4 J. KAMPEN 

its Greeks is set up based on the regularity results of the previous Section 
and the nonlinear integral equation characterizing the free boundary func- 
tion of Section 4. The description is distributed over four subsections. In the 
first subsection subproblenis related to the nonlinear integral equation char- 
acterizing the free boundary function are identified. The second subsection 
recalls some recent results of the WKB expansion of the transition density 
(fundamental solution of linear parabolic equations). The third subsection 
describes the algorithm and the fourth subsection provides a convergence 
and qualitative error estimate analysis. A conclusion is given in Section 7 
where we also indicate future research. Finally, in the appendix we add a 
technical proof needed in Section 5. 

2. General framework. In this Section we shall set up the general 
framework and the technical conditions which are needed for the strong so- 
lutions of the systems of stochastic differential equations which model the un- 
derlyings (and possibly background processes) and the associated variational 
inequalities which determine the value of the American derivative written 
on that underlyings. Let O be an open domain in M" and let T G (0, cxo) be 
the time horizon. We write Z) = (0, T) x O for the domain where O has to 
be specified case by case. For each starting point (t, x) E [0, T] x O, consider 
the following stochastic differential equation 

X*'" = X e O, 
(2.1) 

for continuous drift functions 

(2.2) ^fc : [0,T] xR"^M, k = l,---,n, 
and local volatility functions 

(2.3) aij : [0, T] X M" ^ M i = 1, • • • , n, 

with an M™-valued Brownian motion W = (W^)j=i^... ^m and where X*'^ = 
(X*'^'''^)fc=i ... „. We assume that the functions /i^ and ajk are locally Lipschitz- 
continuous in x, uniformly in t, i.e. for each compact subset K C M" there 
is a constant c (dependent on K) such that 

(2.4) \fik{t,x) - fik{t,y)\ <c\x -y\ 
and 

(2.5) \aij{t,x) - aij{t,y)\ < c\x - y\ 

imsart-aap ver. 2007/12/10 file: aapamerican.tex date: July 7, 2008 



FREE BOUNDARIES AND THEIR GREEKS 5 

for all t G [0,r] x,y £ 0,. The latter assumption implies that (2.1) has a 
unique strong solution for any given filtered probability space ( Q,, J-, (•^t)tg[o,T] ) P 
and Brownian motion 1^ up to a possibly finite random explosion time. 
Therefore we add the assumption that 

(2.6) pI sup |X*'^| <oo) =1. 

\se[t,T] J 

Then theorem V. 38 of [2G] implies that (2.1) has a strong solution. In 
this context we consider a stochastic volatility market model system with a 
process St modeling the price of n assets and an d-dimensional background 
process Y driving the local volatility. Then for each starting point (t, x) £ 
[0, T] X O C M!p X M"* the stochastic differential equation 

(5l'-'^y/'^) = (x,y)GO, 

(2.7) ^ = ris,Xi'^)ds + E]=i ak,is,Xl'-,Y^''-)dWi, 

dY'.'y'' = ui{s, Y,''^)ds + EU Pijis, Y,''-)dWi 

has a unique strong solution, if the drift functions 

r : [0, T]xW^ ^ M, 
(2.8) 

I// : [0, T] X M'^ ^ M, l = l,--- ,d, 

and the (volatility of) volatility functions 

(2.9) A, :[0,r] xM'^^M, l£{l,...,p}, je{l,...,4 

satisfy conditions of form (2.5) and (2.6) respectively. 

Then the price V{t, x, y) of an American option with payoff (p is given by 



(2.10) V{t,x,y):= sup Eg 

reStopjj J,] 



exp 



r(.,5,M</,(r,5*'-'^) 



where Stopr^-pi is the set of all JT^ stopping times with value in [t,T]. The 
measure Q is some equivalent measure to be chosen if the market is incom- 
plete and uniquely determined in a complete market. If the asset process 
S does not depend on a background process Y, then the market process is 
complete, especially S*'^'^ = 5*'^'. For simplicity of notation we deal with 
this situation in the following, noting that everything can be generalized 
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6 J. KAMPEN 

to the more general incomplete market model (2.7) without problems (ex- 
cept the choice of the equivalent measure which is a problem of its own). 
Standard stochastic control theory can now be used to show that the value 
function (t, x) — > V{t, x) satisfies the nonlinear Cauchy problem 

f max {^ + Lu, (p - u\ =0, in [0,r) x M" 
(2.11) 

[ u{T,x) = (l){T,x), in {T} x M" 

in the viscosity sense, where Vij = {era )ij and where 

(2.12) Lu^-Y: %S<S,;j^jg^ + . (Y. S.gs^ - «j 

Global solutions for the latter type of equations in the viscosity sense were 
studied in [1] even in a far more general context. It is clear that (2.11) 
formally implies that 



(2.13) 



(f +Ln)(</.-7z) = 0, 

- n > 0, in [0, T) x W a.e. 

u{T,x) = (f>{T,x) on {T} x M", 

and where global solutions can be obtained under rather mild conditions by 
variational methods (cf. [12] ). The equivalence of this system of variational 
inequalities with optimal stopping problems under relatively mild conditions 
have also be studied extensively (cf. [9] ). The set points (t, x) where the value 
function u equals the obstacle (j) is called the exercise region. The complement 
of the exercise region in the domain D is called the continuation reason. In 
[8] it is proved in a rather general context that the optimal stopping time 
is the first time where the value process hits the boundary of the exercise 
region. Therefore the free boundary surface (or the boundary of the exercise 
region) is crucial information in dealing with American Option. Especially 
in the context of hedging information concerning the Greeks near the free 
boundary is crucial. In the present article we set up a method to obtain the 
Greeks in a general class of multivariate models. Information of the value 
function in the continuation region and its Greeks are then a by-product of 
our method. 
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FREE BOUNDARIES AND THEIR GREEKS 7 

Strict ellipticity of the diffusion matric function (vij) and Holder conti- 
nuity of the functions Vij and r together with bounded Lipschitz continuous 
data (j) is enough to ensure existence of global solutions in viscosity sense for 
(2.11) (cf. [1] for global existence for more general conditions and classes of 
problems), or, with a bit more regularity assumptions (derivatives in L^), 
global existence for (2.13) can be ensured by classical variational methods 
in Sobolev spaces. In the following we shall assume that at least the for- 
mer conditions hold. However in the design of the scheme we shall essen- 
tially assume that coefficient functions are C°° with exponential bounds 
of the derivatives. More precisely we shall assume that for two positive 
constants and all (t, x) and all real n-dimesnional vectors .^ 7^ we have 
< A < Vij{t,x)(,i^j < A and that Vij and r are Holder continuous (with 
exponent a G (0, 1) and bounded. This ensures the global existence of the 
fundamental solution (transition density). 

3. Frontfixing w^ith special focus on American Basket Put op- 
tions. We use a transformation involving the free boundary which is oth- 
erwise known as a frontfixing method. Our method works for a large class 
of options, e.g. for 

• the minimum put with payoff {K — minimi, ^2}) , 

• the spread with payoff {K — {Si — S2)) , 

• the put on an index or index spreads {K — J27=i oaSi) , 

• similar call options in standard models and other standard options. 

Since suitable front-fixing depends on the pay-off we shall restrict the anal- 
ysis to the most popular example the American Put on an index. It will be 
clear from the argument that the other options can be treated similarly. We 
start with the operator 

where Vij = (cr(T-^)jj and r may depend on time t and spatial variables S. 
Let £ denote the exercise region and for each t £ [0, T] let £t denote the 
t-section of the exercise region, i.e. £t '■= {x\{t,x) £ £}. In general beside 
basic standard assumptions on stochastic volatility models introduced in 
the previous Section we shall assume that 

(GG) for each t £ [0, T] G £t and that for a fixed "angle" at 5 = (^i, • • • , 5^), 

i.e. at 

, f S2 Sn 

9s ■-- 




En Q ^ Y^n Q 
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8 J. KAMPEN 

we have one intersection point of the section E-t and the intersection 
point of the ray through which is determined by the angle (^5. 

Remark 3.1. We call (GG) the global graph condition. It is clear that 
this condition is satisfied for the multivariate Black-Scholes model and for 
a considerable class of stochastic volatility models. It is also satisfied by 
market models based on Levy processes. However, the latter type of models 
will be considered elsewhere. 

Hence, the free boundary can be written in terms of the angles in form 
(3.2) {t,ct)s) ^ F{t,^s). 

We consider the transformation 

V' : (0, T)xWl^ (0, T) X [1, 00) x (0, 1)""^ , 
(3.3) 



^(t,5i,--- ,5„) 



TLiEi 



S2 



■Jn 



^ 't::^''''" '^zj-^ 



Note that the spatial part of tl^{{0,T) x M") is homeomorph to the half 
space H>i = {x € M^jxi > 1}. In the following the domain D is the interior 
of the image of V, i-e. D := (0, T) x (1, 00) x (0, 1)""\ We have 



(3.4) 
We get 



(3.5) 



Si = xiF 1 — ^ Xj , Sj = XjXiF. 
i>2 



Ur 



EL^.du_ I 1 



,F d'^u 



F du 



du 



T — pXi Q^^ -I- 2 l^ij "jj Oxidxj ^ l^j "j dxj ^ ' -^1 dxi 



(SCI) u(0, 00, X2-, • • • , Xn) = on xi = 00 



{BC2) u.x,{t,l,X2 



1 ) -^nj 



u{t, 1,X2 



1 ) -^nj 



-K 



on xi = 1 
{BC2,) F{t, X2,- ■■ ,Xn) = K - u{t, 1, a;2, • • • , Xn) 

(IC) u{0,x) =max{K -xi,0} 

Remark 3.2. We include (BCl) as an implicit boundary condition in or- 
der to indicate that (3.5) is equivalent to an initial-boundary value problem 
of the second type on a finite domain (just by suitable additional transfor- 
mation with respect to the variable xi). 
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FREE BOUNDARIES AND THEIR GREEKS 9 

The mixed condition {BC2) follows from the smooth fit condition together 
with {Bd) . 

Remark 3.3. Note that in the context of market models based on Levy 
processes or, more generally, Feller processes the smooth fit condition does 
not hold in general and one has to be careful concerning generalization at 
this point. Especially one cannot expect that a condition of type {BC3) 
holds in the pointwise sense if the smooth fit condition fails. 



In order to det 


iermine the coefficients af, and bf ^ 






Tpdx, _ 5ii-xi 

^ aSi - xi ' 


(3.6) 




a Y- axj a 
aSi ~ ^3 aSi axj ■ 


We observe that 






(3.7) 






It follows that 






E.S, 


a 
aSi 


— l^tj dSi dSi 



0. 



(3.8) 



EC axi a _|_ Y^ / Y^ c 8x1 
i "^^ dSi dxi "•" ^j>2 \2^i "^t aSi 

EC ax I a 
i "^^ dSi dxi 

^IfeT ~ (Ei'S'i) (Ei>2 [dSl] 'F'^^) feT 



^1 axi ■ 



Hence, we have 



i^''^) 



(3-9) r\y:S~]=rx,^^. 
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10 J. KAMPEN 

It is clear that 

(3.10) <=i:-'^iSkSig^^Qs-^, 



d' 






and that 

(3.11) bf = Y:vMS,S,^^^^^. 

In order to determine the latter coefficient function we compute 

^''•^^^ dS,~F XI '-^-^ 

and 

(3.13) ||=i,(i_|:„,_.,)S 

Next, for j > 2 we have 

(3.14) ^x^^yd{^) 



SSj^^fc V ^3;, 



and 



(3 15) ^{^) _ S,i ^ {xj - 6ij){6iiF + xiFijl - 6u) 



dxi xiF {xiFy 



Finally, 



^'■''^ dMs-r^d^^F- ^}'- - ^^^F^j 55^' 

where 
(3.17) 

— —p^U ~ Oil) - l^j>2 F3" 

Here 5ij is always the Kronecker Delta, Fj is short for g^, Fji is short for 

g^ Q^ ■ Now we have determined the explicit form of (3.5). The next step is 
to construct a representation of the solution of (3.5) in terms of convolutions 
with the transition density, i.e the fundamental solution related to (3.5). 
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FREE BOUNDARIES AND THEIR GREEKS 11 

4. Derivation of a system of integral equations characterizing 
the free boundary. For a fixed free boundary function F we may consider 
(3.5) as a standard linear initial value boundary problem of the second type 
and of the form (we write again L> = O x (0, T). 



^'■'^ m + 2^^ ""-d^^ + ^ '^ 9? = ^' '^ ^' 

(4.2) u{0,x) =g{x) onO 

d 

(4.3) —u{t, x) + au{t, x) = h on H^ 

where H := < {t,x) G D\xi = 1 > is part of the boundary of D {H is for 
'hyperplane'). In the standard formulation of the initial value boundary 
problem v denotes the inward normal and Ojj, bi, f, a, h are functions which 
may depend on time t and the spatial variables x. In our case the derivative 
with respect to the inward normal reduces to the partial derivative -^r, 
the function a is the constant function a = —1, h is the constant function 
h = —K, and / = 0. Furthermore, the initial condition for the basket put is 
just g{x) = max{i^ — xi, 0} in the transformed coordinates. Hence equation 
(4.1) simplifies to 



9n 1 " d'^u J^du 

('•') m^2^, ""-d^M + ^ '^5? = °' ^" ^' 

i,j=i 1=1 

(4.5) u{0,x) = inax{K — xi,0} on O, 

d 

(4.6) t;— -"(*, x) - u(t, x) = -K on H. 

OXl 

We represent the solution in terms of convolutions of the fundamental 
solution and an integral equation of Volterra type related to the mixed type 
boundary condition. First we define the fundmental solution p = p{t, x; r, y) 
for t > r to solve for each (r, y) 

«j=l 1=1 

(4.8) u(0,x) =5j^(x) at t = 0. 

The general ansatz for (4.1)-(4.3) then is (recall that xi = {x2, • • • , Xn),yi = 
(2/2, • • • , Un), and dHy = dy2dys ■ ■ ■ dyn), writing H = Hq x (0, T) and D = 

imsart-aap ver. 2007/12/10 file: aapamerican.tex date: July 7, 2008 



12 J. KAMPEN 

O X (0, T) is 

n(t,x) = ^^^H^p{t,l,xi;T,l,yi)(t){T,l,yi)dHydT 
(4.9) 

+ /o ?*(*' ^; 0' y)9{y)dy - /o /o p(*' ^; '^' y)/(T, y)dydT, 

where the boundary condition -§^u{t, x)+a{t, x)u{t, x) = h{t, x) then reduces 
to the Volterra type equation 

l^{t,x) =T{t,x)+ 

(4.10) /o Jh^, im^Pit^ 1> xi,T,yi) + a{t, 1, xi)p{t, 1, £i; r, 1, yi)) 

X(j){T,l,yi)dHydT, 
with 

r(i, x) = Jq ^p{t, x; 0, y)g{y)dy 

(4.11) +a{t,x)Sop{t,x;Q,y)g{y)dy 

-a{t, x) Jq Jq p{t, x; r, y)/(r, y)dydT 

—h{t,x). 

With a = —1, h = —K, and / = 0, and g{x) = maxjii' — a;i,0} this 
simphfies to 



(4.12) 

where 

(4.13) 

with 

(4.14) 



u{t,x) = jQj^^p{t,l,xi;T,l,yi)cl){T,yi)dHydT 
+ lo P(*' ^; 0' y) ™ax {K - yi, 0} dy, 



^^{t,x) =T{t,x) + 

/o /i^o (air^(*' ^' ^1' '^' ^i) + ^(*' ^' ^i; '^' 1' yi))'?^(^' 1' yi)dHydT, 



T{t,x)= Jq -^p{t, x;0,y) max {K -yi,0}dy 

- IoP(^^ ^; 0' y) max{/^ - yi, 0} dy + iT. 
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FREE BOUNDARIES AND THEIR GREEKS 13 

We apply this to initial value boundary problem (3.5) in order to obtain a 
system of integral equations for the free boundary. First, from the boundary 
condition (BC3) in (3.5) 

(4.15) Fit, X2,--- , Xn) = K - u{t, 1, 2:2, • • • , Xn) 

we get the integral equation 

F{t,X2,--- ,Xn) = K -u{t,l,xi) = 

(4.16) K - Jq J^^ p{t, 1, f 1; T, 1, yi)0(T, 1, yi)dHydT 

+ /oP(*'2^;0'y)max{i(:-yi,0}(iy. 

Note that the function (j) which is itself solution of a linear integral equation 
and can itself be represented in terms of the fundamental solution. Define 

1 d 

(4.17) i^Zpit, x; T, y) = Q^PPit, x; 0, y) - prit, x; 0, y), 

and for r > 1 define 

(4.18) Z''p^^{t,x,T,y) = / Zp{t,x,a,l,zi)Zp{a,zi,T,y)dHyda 

Jo J Ho 

Then 

-^ °° ft f 

(4.19) -cP{t,x)=r{t,x) + J2 / ZUt,x,T,l,yi)dHydT 
^ ~iJo J Ho 

We summarize 

Theorem 4.1. Assume that (GG) and the assumptions of the general 
framework hold. Then the free boundary surface function F of the multivari- 
ate American basket put option (with weights normalized to 1 w.l.o.g.) in the 
transformed coordinates 

^l) : (0, T)xWl^ (0, T) X [1, oo) x (0, 1)""^ 

(4-20) ^„^^^ ^^ 

V'(^; Si,- " ) Sn) — y^ — ^'^ — 1 y^" ^ g_ ) ■ ■ ■ ) y^'i" g. 1 

is solution of the integral equation 

F{t, X2,--- ,Xn) = K - u{t, 1, f i) = 

(4.21) K - /o* f^o PF{t, 1, xi; T, y)^{T, y)dHydT 

+ loPFii^ 1, ^i; 0, y) [K - yi)+ dy, 
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14 J. KAMPEN 

where 

1 °° ft f 

(4.22) -^{t,x)=T{t,x) + Y, / Z'^{t,x,T,l,yi)dHydT 

2 ^./O J Ho 

with 

1 d 

(4.23) 2'^^(*' ^' '^' y^ = 'dx~^^^^' ^' ^' ^'' ~ -^■^*-*' ^' ^' ^''' 

(4.24) Z]+\t,x,T,y)= [ f ZUt,x,a,z)Z'^{a,z,T,l,yi)dHyda, 

Jo J Ho 

for r > 1, and pp is the fundamental solution of 

/, „^N Ft du 1 v^ F f^^"" sr^ipdu du 



«J " ^ 3 ^ 



Remark 4.2. Note that the fundamental solution of (4.25) can be given 
in the form of a Levy expansion, i.e. 

(4.26) p^{t,x,T,y) = Z^{t,x,T,y)+J J Z^{t,x,s,z)<^^ {s,z,T,y)dzds 
where 

oo 

(4.27) $^(t, X, T,y) = Y^ (L^ZUt, x, t, y) 

n=l 

(4.28) 

{L^Z)i{t, x, r, y) = L^Z^{t, x, r, y), 

{L^Z)n+i{t, X, r, y) = J^ f^ L^Z^{t, x, s, z){L^Z)n{s, z, r, y)dzds, 



where Z^{t,x,r,y) = -^^=L=^^detA~/ir,y) exp (^_ "'.-(-.^)(^-^_^;)fe-^.) 

with detAp {T,y) = a^p{T,y) being the inverse of {afAr^y)). The Levy ex- 
pansion is not an efficient tool in order to compute the fundamental solution, 
but we shall use it for proving existence and uniqueness of the free boundary 
function below. 
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5. Global existence and uniqueness of solutions. From the rela- 
tion 



(5.1) 



F(t,xi) = K-u{t,l,xi) 



(recall that (x2, • • • , x„) = xi) it is clear that uniqueness of u implies unique- 
ness of F and vice versa. Moreover, the initial value (free)-boundary problem 
(3.5) can be written as a nonlinear initial value boundary problem as in (5.2) 
below. Indeed, let 

u {t, xi) = u{t, l,a;i) 

and let us denote by ajt- and 6f the coefficient functions obtained from af^ 
and hj substituting F hy K — v} and the derivatives of F by the derivatives 
of —u}. Then we have 

(BCl) u{t, oo, X2, • • • , Xn) = on Xi = oo. 



(5.2) 



{BC2) U;,^(t,l,X2,--- ,Xn) -n(t, 1,X2, 



(IC) u{0,x)=max{K -xi,0}. 



c„) = -K, 
if xi = 1 



It is clear from our representation that once regularity of the free boundary 
surface is guaranteed, global existence of the solution can be obtained from 
(5.2). We shall prove regularity of the free boundary function by using the 
integral equation in theorem 4.1. The idea is to look at this integral equation 
as a map of a suitable Banach space where the free boundary function is 
obtained as a fixed point. In order to state the theorem we recall the defi- 
nition of Banach spaces related to Schauder estimates. Schauder boundary 
estimates turn out to be appropriate since the problem in (3.5) is equivalent 
to a problem on a bounded domain with Dirichlet conditions on a trans- 
formed hyperplane H^o- Including infinity into the topology we denote the 
boundary of (3.5) by 



(5.3) 



dD = HiUHooUlUFUS^. 



ang 



and, for convenience, we define 

(5.4) S:=H^UH^U 5ang 
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with 

Hi:=H= {(t,x) G (0,r) xO\xi = l} , 

Hoc ■■= {{t,x) G (0,r) X 0\xi = 00} , 

(5.5) / = |(t, x) G [0, T] X 0\xi > 1 and t = 0} , 

F = {(i,x) G [0,T] X 0\xi > 1 and t = t} and 

Sang = {{t, x) G [0, T] x 0|xi G {0, 1} and i > 2} . 

We assume that (3.5) has been transformed to a finite domain (a suitable ad- 
ditional transformation with respect to the variable xi s easily constructed) 
in the sense that Hoo is transformed to H^ with d > 1. Therefore, in the fol- 
lowing we shall be able to use some results on bounded domains Df, C M"+^. 
Abbreviating £- =: D^^, g^ =: L»^,^^-, | =: A- 

For any points P G S" and Q (^ Df, define the distance from Q to 5 by 

(5.6) dQ:= inf d(P,Q), 

where d{P, Q) denotes the Euclidean distance between P = {t,x) and Q = 
{s,y), i.e. 



(5.7) d(P,Q) = ^|x-y|2 + |t-s|2. 

Furthermore, for any two points P,Q G Db, integers m, > 1 and any function 
V : Db ^f ^ let dpQ = min {dpjdq} and define 

(5.8) 6S^{d-f) := lubp,QeD,rf^^° '"^^^"J^y^ 

and 

(5.9) \v\a' = sup \v{t,x)\+6^^{v) 

Let 



n 



(5.10) big, = bi^" + Y. \dD,Ma^ + J2 M'^x^x.^la" + |d'A^lf^ 
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and define 

(5.11) C2+a{Db) ■.= {f:D^ nf\2ta < °o} 

Recall that 

Proposition 5.1. If D^ C M" is bounded, then C2+a{Db) is a Banach 
space. 

Next we recall some classical results on a priori estimates by Schauder. 
Consider the first initial-boundary value problem on a hypercube D: 

{ ^-\l:^,a^,{t,x)^^-Y.^h{t,x)^^+c{t,x)w = f{t,x) 

(5.12) J inDfcUF, 

y w{t, x) = h{t, x) onlUH X (0, T) U H^ x (0, T), 

where H = <x ^ Dh\xi = 1 r , H^ = Ix €^ Dh\xi = d> (for some d > 1, and 
/ and F are the boundaries of D at t = and t = T respectively. Assume 
that the following assumptions are satisfied: 

(A) the coefficients {t,x) — > aij{t,x) satisfy an ellipticity condition, i.e. 
there exists a constant C > such that for any (t, x) G Dh 

(5.13) J^aij{t,x)C^^,>c\C\^. 

(B) The coefficient functions (t,x) -^ aij{t,x) and {t,x) — > bi{t,x) and 
{t,x) — > c{t,x) are locally Holder continuous (exponent a), i.e. there 
exists a constant C > such that 

(5.14) \aij\a < C, \dbi\a < C, {d^cla < C 

(C) the function / is locally Holder continuous (exponent a), i.e. there 
exists a constant C > such that 

(5.15) \d^fU < C 

Remark 5.2. We say that (A), and (B) of the present section hold for 
the domain D for the problem (3.5) if (A) and (B) hold for the equivalent 
problem defined on a finite domain Df, where H^ has been transformed to 
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Definition 5.3. We say that S has the outside strong sphere property 
if for every Q = {t,x) £ S there exists a ball B with center {tm,Xm) such 
that BnDi) = Q and |x — x^l > ^{Q) > for all {t,x) £ D}, and |t — tm| ^ £■ 

We observe 

Proposition 5.4. The boundary surface H has the outside strong sphere 
property. 

Proof. Trivial. 

Now a classical result on Schauder estimates states that 

Theorem 5.5. Assume that (A), (B), and (C) hold and assume that S 
has the outside strong sphere property. Then for any continuous function h 
on I U H there exists a unique solution of the (5.12), and u G C2+a{Db)- 

Now we can establish the main result. 

Theorem 5.6. Assume that the assumption of the general framework, 
the assumptions (A), and (B) of the present section hold for the domain 
D (cf. remark 5.2) , and the global graph condition (GG) holds. Then free 
boundary function F ■]Q,T[x]Q,l[''-^^M. in (3.5) isinC2+a{D). 

Proof. The main observation here is that we can use a priori estimates of 
the first initial-boundary value problem to study the regularity of a class of 
linear second initial-boundary value problems. We first show that there is a 
function 

G : C2+a{D) ^ C2+a{D) 

(5.16) G{v){t,x) = K - jQjjjpy{t,l,xi;T,l,yi)(py{T,l,yi)dHydT 

+ loPvit A, xi;0, y) {K -yi)'^ dy, 

i.e. we show that G as specified in (5.16) is indeed a function from C2+a{D) 
to C2+a{D). Note that the classical theory of the fundamental solution (Levy 
expansion) immediately implies that {t,x) — > G{v){t,x) is continuously dif- 
ferentiable with respect to time and twice continuously differentiable with 
respect to the spacial variables for v G C2+a{B>). However, both integrals in 
(5.16) have an interpretation as an inititial-boundary problem of type (5.12) 
(note the implicit data on the boundary) and application of theorem 5.4. 
leads to the conclusion that G is indeed a map of C2+a{D). 
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Lemma 5.7. For all t G [0,r], x e D and ■wi,W2 G C2+aiD) there is a 
constant L < oo such that 

(5.17) \G{wi) - G{W2)\^^^ < Lt\wi - W2\^^^. 

Proof. The proof is not difficult but technically cumbersome. We indicate 
the basic idea. The proof is based on the Levy expansion of the fundamental 
solution. For fixed Wi,i G {1,2} consider the Levy expansions (cf. remark 
4.2.) of the fundamental solution 



(5.18) p''^{t,x,T,y)=Z'"^{t,x,T,y)+l J Z^{t,x,s,z)^''^{s,z,T,y)dzds 



ft 
Id 
where 



(5.19) $"''(t,x,r,y) = lim ^Z'it,x,T,y) := lim ^(L^Z)„(t, x,r, y), 

n=l 

and its derivatives up to second order. Note that the second spatial deriva- 
tives have the form 

p'''it,x,T,y)= ^Z-''it,x,r,y) 
(5.20) 

+ IrlDdS^^^it,X,S,z)^^'{s,Z,T,y)dzds, 



the first time derivative has the form 

d_ 
dt 



p^^{t,x,T,y) = §rZ^^{t,x,T,y) - <^^'{t,x,z,T,y) 



(5.21) 

+ Ir Id m^^it, ^, s, z)<^^^ {s, z, T, y)dzds, 



and the first order spatial derivatives are of the form 

d 

dxi' 



p^^{t,x,T,y)= J^Z-'(t,x,r,y) 



(5.22) 

I r* r 

dxi 



+ Ir Id m^Z^ii^ ^> ^' ^)^"' (^. ^> r, y)dzds. 



Inserting these expressions into \G{wi) — G{w2)\2a ^-^d using triangle in- 
equalities and induction over m leads to the desired result. D 

We conclude that the function G : C2+a{D) — > C2+a{D) is a contraction 
and its unique fixed point, and the back transformation with respect to xi of 
this fixed point is the free boundary function F, which is hence in G2+a{F)). 
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6. Iterative schemes for the free boundary surface and the Greeks. 

The contraction map G in the regularity proof above leads to various versions 
of iterative schemes for the free boundary surface F, its time derivative and 
its spacial derivatives up to second order. We describe an algorithmic scheme 
which has several possible realizations (probabilistic and PDE-schemes) of 
its subproblems. We shall describe the probabilistic scheme in more detail. 
PDE-schemes and issues of implementation will be considered elsewhere. 

6.1. Splitting scheme (splitting the mixed boundary problem at each itera- 
tion step). In the proof of theorem (5.6) we observed that the free boundary 
function F is a fixed point of the intergal equation 

F{t,xi) = K - J^JjjPFit,l,xi;T,l,yi)(t)F{T,l,yi)dHydT 
(6.1) 

+ loPPii, x'l 0, y) max {K - yi, 0} dy. 

along with a linear integral equation for (pp- Now given a function F G 
C2+a{F)) we may first solve for the linear integral equation for (f)F. This is a 
Volterra integral equation well-studied in the literature numerically (cf. [27] 
for a probabilistic treatment). Next we analyze (6.1). First, the term 

(6.2) pF{t,x;0,y)max{K -yi,0}dy 

Jd 

can be interpreted to be the solution of 

(6.3) ; 

(IC) u{0,x)=max{K -xi,0} 
on D with natural boundary conditions on H. The interpretation of 

(6.4) / / pF{t,l,xi;T,l,yi)(l)F{r,l,yi)dHydT 

Jo Jh 

as a solution of a Cauchy problem with initial condition equal to zero and 
a source related to (pF of reduced dimension n — 1 is more theoretical in 
the sense that a statement of the related lower dimensional Cauchy prob- 
lem is in general only possible if the solution u is known. This is a first 
reason why an efficient computation of the transition density (fundamental 
solution) is desirable, and it also indicates that our algorithm is intrinsi- 
cally probabilistic to some extent in the sense that it cannot be completely 
realized without referring to the transition density. A second reason is the 
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following. If we know the fundamental solution pp in (6.24) and in (6.3) (for 
fixed F), then we can obtain the expressions pF,t, PF,Xi,i G {2, • " " i^}i and 
PF,xix ihj £ {2, • • • , ra} for the next iteration step n+1 by differentiation un- 
der the integral, i.e given the approximation of the free boundary surface of 
the ra-th iteration step F" and its derivatives occurring in pp^ we get the ap- 
proximation of the free boundary surface of the n + 1-th iteration step F"'+^ 
and its derivatives Fj"+\ i^"+\i G {2, • • • ,n}, and F^+'^,i,j G {2, ••• ,n} 
via 

^F"+i(t,xi) = -§-^ J^ JjjPF«{t,l,xi;T,l,yi)(l)F"{T,l,yi)dHydT 
+ Id mPP'" (*' ^; 0, y) max {K - yi, 0} dy 



(6.5) 



^F"+i(t, xi) = - /o Xf^ §iVF^{t, 1, xi; T, 1, ?/i)</'F" (t-, l,yi)dHydT 
a 



+ /d ^Pi^" (*' ^; 0' y) ™ax {K - yi, 0} dy 



^F"+i(t,xi) = -/o7^g^Pi.n(t,l,£i;r,l,yi)x 






6Fn(r, l,yi)dHydT + /^ 3^;^^PF^{t, x; 0, y) max {K - yi, 0} <iy 



Note that 

^F"+^(t, xi) = - /o /j:^ ftPF"{t, 1, xi; r, 1, yi)(/'F" (t, l,yi)dHydT 

(6.6) + /j:^ pF"(i, 1, ^1 ; T, l,yi)(pF" {t, l,yi)dHy 

+ Xd WtPF"^ {t,x;0,y) max {/sT - yi , 0} dy , 

giving a formula for ^F""*"^ in terms of the transition density pp"- As we 
observed in [15], approximating derivatives of value functions by derivatives 
of analytical approximations of transition densities in convolution represen- 
tations of the value function is numerically efficient and allows for error 
estimates in strong norms by using a priori estimates of Safonov type (cf. 
[19]), and applying them to analytical representations of the transition den- 
sity and its approximations (cf. [15]). The WKB-expansion of the transition 
density is such an analytical representation. 

6.2. Use of WKB- expansions of the fundamental solution at each iteration 
step. We review some recent research on the fundamental solution (tran- 
sition density), i.e. some results concerning WKB-expansions of parabolic 

imsart-aap ver. 2007/12/10 file: aapamerican. tex date: July 7, 2008 



22 J. KAMPEN 



equations (cf. [15, 17], for more details). The solution (t, x, s, y) -^ p{t, x, s, y) 
of the family of parabolic equations 



du 1 Y^ ri 9 tt \^ h c)tt n 

dt 2 ^i,j "^i dxidx, ^i "« dx, ~ ^ 

(6.7) 

p{0,x,0,y) = 6{x-y), 

on a domain D, parameterized by y G M", and where 5 denotes the Dirac 
delta distribution and the diffusion coefficients aij and the first order coeffi- 
cients bi in (6.7) may depend on time t and the spatial variable x. Without 
loss of generality and for simplicity of notation we consider the case where 
the coefficients depend on the spatial coordinates. Note that the time coor- 
dinate in the coefficients can be treated as an extra spatial coordinate and 
the resulting degenerate parabolic equation belongs to a class of so-called 
projective parabolic equations which are subject to all the following results 
(cf. [18]), In the following let 5t := t — s, and let the functions 

{x, y) -^ d{x, y) > 0, (x, y) -^ Ck{x, y), k>0, 

be defined on [0, T] x M" x [0, T] x M". Then a set of (simphfied) conditions 
sufficient for pointwise valid WKB-representations of the form 

(6.8) p{t, x,s,y) = -^=^ exp I - ^^^ + £ c^ (x, y)6t^ j , 

is given by 

(WKBl) The operator L is uniformly elliptic in M", i.e. the matrix norm of 
{aij{t, x)) is bounded from below by A > and from above by A > A, 
uniformly in x, 

(WKB2) the smooth functions {t,x) — > aij{t,x) and {t,x) -^ bi{t,x) and all 
their derivatives are bounded. 

For more subtle (and partially weaker conditions) we refer to [17]. If we add 
the uniform boundedness condition 

(WKB3) there exists a constant c such that for each multiindex a G N", and 
for all I < i,j,k < n, 

(6.9) ^^^^) ' 9^(^) ^^^^P(^l^l> 

then the Taylor expansions of the functions d and c^ around y E M" are equal 
to d and Ck,k > globally, i.e. we have the power series representations 

Ck{x,y) =J2aCk,a{y)^^'^^ ^ > 0, 
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where 6x := x — y. Note that (C) is impUed by the stronger condition that 
ah derivatives in (6.9) have a uniform bound. Summing up we have the 
fohowing theorem: 

Theorem 6.1. If the hypotheses (A),(B) are satisfied, then the funda- 
mental solution p has the representation 

(6.11) p{6t, X, y) = 1-^ exp "^^^ + E Cfc(a^> y)<5i^ 

where d and c^ are smooth functions, which are unique global solutions of the 
first order differential equations (6. 12), (6. 13), and (6.15) below. Especially, 



n ^ , , ^ , d 



•2 



{6t, X, y) -^ dt liipidt, X, y) = --6t ln{27r6t) - — + Y.'^k{x, y)6t^ 

^ ^ k>0 



is a smooth function which converges to —-^^ cls 5t\^ 0, where d is the Rie- 
mannian distance induced by the line element ds"^ = J^ij^^ij dxidxj, where 
with a slight abuse of notation (a^- ) denotes the matrix inverse of {oij). If 
the hypotheses (A), (B) and (C) are satisfied, then in addition the functions 
d, Ck,k>0 equal their Taylor expansion around y globally. 

The recursion formulas for d and c^, k > are obtained by plugging 
the ansatz (6.8) into the parabohc equation (6.7), and ordering terms with 
respect to the monoms 5t* = (T — i)* for i > —2. By collecting terms of 
order 6t~'^ we obtain 

(6.12) d2 = i^4^a,,4^., 



4 



«j 



where d^ denotes the derivative of the function d^ with respect to the 
variable Xk, with the boundary condition d{x,y) = for x = y. Collecting 
terms of order St~^ yields 



(6.13) - ^ + hd' + ^ E ( E K(^) + «.-^(^)) ^ 

where the boundary condition 

(6.14) co{y,y) = --\n^det{aij{y)) 



S("'^)=°' 
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determines cq uniquely for each y € M". Finally, for A; + 1 > 1 we obtain 

(6.15) 

with boundary conditions 

(6.16) Cfc+i(2;,y) = i?fc(y,y) if x = y, 

i?fc being the right side of (6.15). In [18] it is shown how the function cP can 
be approximated in regular norms if only (WKBl) and (WKB2) are satisfied. 
The technique can be combined with the regular polynomial interpolations 
developed in [10]. 

6.3. Description of algorithm (including to sparse grids and weighted Monte- 
Carlo versions). We can solve (3.5), (or, similarly) (5.2) by an iterative 
numerical procedure. 

For simplicity of description we assume that the Lipschitz constant k = 
k{Dt) of G above with respect to \.\2a is less than 1. This is true in any 
case for t small enough. It is clear that iteration in time of the scheme below 
leads to a global time scheme. 

(Stepl) Solve 



1 _ 1 Y^ F d^^i I v^ y,F duj , dul 

^l,r - 2 ^ij ^ijdxidxj + ^j "j dxj + "^^^ dxi ' 

(6.17) ; 

(IC) u}(0,x) =max{A'-xi,0}, 
where ajj '•= o-fj ^] •= ^f for F = 1, i.e. compute 

(6.18) u\{t,x):= I pi{t,x;0,y)msix{K -yi,0}dy, 

Jo 

where pi(t, x; 0, y) := prit, x; 0, y) with F = 1. 
Next solve for (pi in 

i(^i(i,x) = T{t,x)+ 
(6-19) Jq Jfj^^{^pi{t,l,xi,T,yi) +pi{t,l,xi;T,l,yi))x 

(j^{T,l,yi)dHydT, 
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where c/)^ = (p for F = 1. 
Next compute 

(6.20) u\{t,xi):=l / pi{t,l,xi]T,l,yi)(l)F{r,l,yi)dHydT 

Jo JHo 

and 

(6.21) Fi{t,xi) = K - ul{t,xi) + u\{t,x). 

Next compute the time (first order) and spatial derivatives (up to 
second order) of Fi according to the formulas (6.5) and (6.6) for F = 1. 

Remark 6.2. The use of WKB-expansions of pi seems the most ef- 
ficient method. In [15] it was observed that a WKB approximation of 
the transition density including the term ci allows pricing of options 
and its derivatives with respect to the underlyings with maturity of 
10 years in LIBOR models of dimension n = 20 with one time step 
(beating all concurrent methods). Moreover, efficiency and accuracy is 
kept for the derivatives because of explicit WKB-approximations. 

(Step2) Having computed F„, For n > 1 compute Fn+i as follows. 
Compute 

(6.22) <(i,x):= [ pF^{t,x;0,y)max{K-yi,0}dy. 

Jo 

Next solve for 6"' in 



i(/)"(i,x) = r(t,x)+ 
(6-23) /o iHoim^PF^ {t, 1, xi,T,yi)+ pF„ {t, 1, xi;t, 1, yi)) x 

4>'^{T,l,yi)dHydT, 
where (j/^ = cp ". 

Next compute 

(6.24) u^{t,xi):= f f pFAtA,^i;r,l,yi)r{r,l,yi)dHydT 
Jo JHo 
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and 

(6.25) Fn{t, xi) = K- ul^{t, xi) + < (i, x). 

Next compute the time (first order) and spatial derivatives (up to 
secons order) of -Fn+i according to the formulas (6.5) and (6.6) for 

F = Fn. 

(Step3) The functions F„ form a Cauchy sequence in the Banach space C2^^- 
Hence we iterate the step 2 until 

(6.26) \Fn+m - F„,|i^„ <Y^\F,- 1|^, < e 

for a prescribed e > 0. 

Remark 6.3. The initial step 1 and the iterated step 2 performed with 
WKB approximations may be done with different methods based on the 
dimension of the problem: 

Ml For problems up to dimension 3 computation grids such as UG may 
be used. 

M2 For higher dimensional models up to dimension n = 5, 6 sparse grid 
techniques may be used (cf. the analysis in [14] and [24]. At present -at 
least to my knowledge- nobody came up with stable numerical sparse 
grids solutions of higher dimension for linear parabolic problems with 
a comparable complexity. 

M3 In any case Monte-Carlo realizations of the algorithmic scheme pre- 
sented above are possible. Recently weighted Monte-Carlo schemes 
have been developed in [13] and [15], where the improved estimators 
established in the latter article allow also to deal with highly peaked 
densities which occur especially for small time or small volatility. 

Remark 6.4. For the implementation of the algorithm some recent re- 
sults on regular polynomial implementation (cf. [17]) and the computation 
of the Riemannian metric (cf. [IS]) in regular norms are needed. 

Remark 6.5. Let us look at the case of higher dimension and proba- 
bilistic realizations of the algorithmic scheme presented above. Given a free 
boundary approximation function F„ the Monte-Carlo estimators used for 
the computations of (6.22) and of (6.24) are based on the formula 

(f\97\ ^""^ (t\ - f ^" PFSi^x,g{x,i))u{g{x,i)) 

^ ' dx'^^' dx- ^{t,x,g{x,0) 
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and -^^ are spatial derivatives up to secons order or a similar formula for the 
time derivative. Here ^ be an R"-valued random variable on some probability 
space with a density \{z) ^ for all z, and the regular (at least twice 
continuously differentiable) map C,^ := g{x,$,) satisfies \dg{x, z)/dz\ / 0, 
and has density (j){x, •) on R" . The corresponding Monte Carlo estimator is 

9^" M^^^dx cl){x,g{x,mO) 

In [15] it is shown that the latter estimator is of bounded variance even for 
small time or volatility. Furthermore the extension to Bermudean options 
in [15] can be transferred in order to compute the Greeks for the equation 
(6.22). Moreover, Monte Carlo methods for computing (6.23) and (6.23) can 
be found in[27] 

6.4. Convergence and error analysis. We confine ourselves to a concep- 
tual error analysis for the higher dimensional case where the Monte-Carlo- 
estimators in (6.28) is used. We content ourselves with an analysis which 
shows the convergence of the algorithm described above, i.e. the analysis is 
qualitative in the sense that the errors related to the projections of analyt- 
ical approximation functions (such as WKB-expansion approximations) on 
grids (or sparse grids) and the numerical error induced by the computation 
of the Volterra integral equation are neglected. A more detailed analysis is 
beyond the scope of the present work and will be provided elsewhere. 

The algorithm described above describes an iteration at each time step 
which involves the solution of an integral equation and two problems (6.22) 
and (6.24) both of which can be interpreted to be of parabolic first initial- 
boundary value problem type. For the computation of these problems we get 
a truncation error of the WKB-expansion and a Monte-Carlo error induced 
by replacing integrals by finite sums. In order to see how this turns out let 

(6.29) ^^^^^(t,.) =i?£;p^^^.'(t,x,C)^, and 
dx°'\''^^) — M dx°'P '-'''•^'™ ''/(/>(*,„()' 

where p^^^^^ denotes the /th order WKB- approximation 

(6.30) Pi{t,x,T,y) = -^=i=^exp (-^-^ + J2ck{x,y)6A . 
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Then the principal error in each iteration step is estimated by 



3S^v*'^) 



d°'I(f ^^ 



< 



^^{t-,x) 



dl^ 



-it,x) 



(6.31) 



+ 



Qa:jWKB,l _ Q^ 

ax°' y-^) dx" v*-' -^ 



leading to the error for the computation of (6.22) and (6.24). One can show 
that the convolution ui of (6.30) with bounded continuous data is integrable 
for / large enough by comparison with the classical estimates for the Levy 
expansion of the fundamental solution. For large I one gets 



(6.32) 



\u{t, X) - Ul{t, x)|i+5/2,2+<5 G 0(t' 2 ). 



using Safonov estimates cited in [19], and where u is the true value function, 
i.e. the convolution with the transition density p (cf. [15] for more details). 
In each iteration step then we compute functions F^ and its derivatives -F", 
F^.,i £ {2, • • • ,n}, and F'^.^i^j G {2, • • • ,n} on a grid which determines 
the coefficient functions 



pn 

a '-J 



and h i 



on the same grid. The approximating Monte-Carlo sums then approximate 
the coefficients with respect to the expectation Holder norm (cf. [15]). Next 
we can invoke stability results of parabolic equations with respect to coef- 
ficients and with respect to initial data for continuous affine interpolations 
of a 'J and b J using the approximations on the grid points. This shows the 
convergence of the probabilistic scheme. 



7. Epilog. We have proved regularity of the free boundary surface for 
a considerable class of relevant market models, and we have set up in detail 
a scheme for the important and difficult problem of computing the Greeks 
for American type options even in the context of higher dimension. On 
the way we have constructed a nonlinear integral equation which charac- 
terizes the free boundary in the multivariate case. For the implementation 
of the scheme in more general situation an implementation of the compu- 
tation of the Varadhan metric in regular norms and of the drift functions 
is needed. This is possible using regular polynomial interpolation (cf. [16] 
and the analysis of the eikonal equation characterizing the Varadhan metric 
global approximation of its solution (cf. [18]). The scheme established is very 
flexible. Note that PDF-scheme realizations are possible beside Monte-carlo 
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realizations especially for lower-dimensional models. Detailed error analysis 
for different realizations of the scheme and their implementations is certainly 
of interest. Furthermore, extension to jump-diffusion models is possible and 
will be considered elsewhere. 

References. 

[1] Amadori, a. L., Nonlinear Integra-Differential evolution problems arising m option 
pricing: a viscosity solution approach, Differential and Integral equations, (2003), Vol. 
16, no. 7, p. 787-811. 
[2] Belomestny, D., Milstein, G.N., Monte-Carlo evaluation of American options 
using consumption processes. International Journal of theoretical and applied finance, 
Vol. 9, Nr. 4 , 455-481, 2006. 
[3] Blanchet, a., Dolbeaut, J., Monneau, R., On the one-dimensional parabolic ob- 
stacle problem with variable coefficients. Progress in Nonlinear Differential equations 
and their Applications, Vol. 63 Birkhauser Basel, p. 59-66, 2006. 
[4] Broadie, M., Detemple, J., The valuation of American Options on multiple assets., 

Math. Finance 7, (1997), no. 3 S. 241-286. 
[5] Caffarelli, L., The regularity of free boundaries in higher dimension, Acta. Math. 

139 (1977), no. 3-4 S. 155-184. 
[6] Caffarelli, L., Petrosyan, A., and Shangholian, H., Regularity of the free 
boundaries in parabolic potential theory, J. Amer. Math. Soc. 17 (2004), pp. 827-869. 
[7] Elie, R., Fermanian, J.D... and Touzi, N., Kernel estimation of Greek weights 
by parameter randomization, Annals of Applied Probability 2007, Vol. 17, No. 4, 
1399-1423 
[8] El Karoui, N., Les aspects probabiliste du controle stochastique. Ninth Saint Flour 
Probability Summer School, (1979), Vol. 876, Lecture Notes in Mathematics, 73-238, 
springer, Berlin. 
[9] Fleming, W. H., Soner, H.M., Controlled Markov processes and viscosity solutions., 
Springer, 2006. 
[10] Friedman, A., Parabolic variational inequalities in one space dimension and smooth- 
ness of the free boundary., J. of Functional Analysis 18, (1975), no. 3 S. 151-176. 
[11] Friedman, A., Partialdifferential equations of parabolic type, Prent. Hall, (1964). 
[12] Friedman, A., Variational principles and free boundary problems, Robert E.Krieger 

Publishing Co. Inc., Malabar, FL. 
[13] Fries, C, Kampen, J.: Proxy Simulation Schemes for generic robust Monte Carlo 
sensitivities, process oriented importance sampling and high accuracy drift approxi- 
mation (with applications to the LIBOR market model). Journal of Computational 
Finance, Vol. 10, Nr. 2, 97-128, 2007. 
[14] Griebel, M., Oeltz, D., A sparse grid space-time discretization scheme for parabolic 

problems. Computing 81, (2007), S. 1-34. 
[15] Kampen, J., Kolodko, A., Schoenmakers, J., Monte Carlo Greeks for financial 

products via approximative transition densities, to appear in SIAM J. Sc. Comp. 
[16] Kampen, J., Regular polynomial interpolation and approximation of global solutions 

of linear partial differential equations, WIAS preprint 1230, 2007. 
[17] Kampen, J.: The WKB-Expansion of the fundamental solution of linear parabolic 
equations and its applications . book manuscript, accepted by Mem. Am. Math. Soc. 
(downloadable at http://papers.ssrn. com/abstract=914122, Soc. Sc. Res. Net), 2006. 
[18] Kampen, J., Computing the length of a geodesic on a Riemannian manifold in regular 
norms, in preparation. 

imsart-aap ver. 2007/12/10 file: aapamerican.tex date: July 7, 2008 



30 J. KAMPEN 

[19] Krylov, N.V., Lectures on Elliptic and Parabolic Equations in Holder Spaces, Grad- 
uate Studies in Mathematics, Vol. 12, American Mathematical Society, 1996. 
[20] McKean, H.P.,J., Appendix: a free bundary problem for the heat equation arising 

fom a problem in mathematical economics, Indust. Manage. Rev., 6: 32-39. 
[21] Milstein, G.N., Reiss, O., Schoenmakers, J. M. A New Monte Carlo methods 

FOR American options Int. J. of Th. Appl. Fin., 7(5), 591-614, 2004. 
[22] Muthuraman, K., a moving boundary approach to American Option Pricing., 

Preprint, ssrn.com/abstract=910882. 
[23] Nielsen, B. F., Skavhaug, O., Tveito, A., Penalty and front-fixing methods for 

the numerical solution of American option problems., The Journal of Computational 

Finance, 5 (4), 69-97. 
[24] Petersdorff, T. v., Schwab, Ch. Numerical solutions of parabolic equations in 

high (iimen5ion,Mathematical Modelling and Numerical Analysis, 2004, vol. 38, no 1, 

p. 93128. 
[25] Petrosyan, a., H. Shahgholian, H. Parabolic Obstacle Problems applied to fi- 
nance: A free-boundary-regularity approach. Preprint 2007. 
[26] Protter, p.. Stochastic integration and Differential Equations, Springer, Berlin 

2004. 
[27] Sabelfeld, K., Monte Carlo methods in boundary value problems. Springer, Berhn 

1991. 
[28] Shahgholian, H. Free boundary regularity close to initial state for parabolic obstacle 

problem, TAMS 2008 (in press). 
[29] VAN Moerbeke, p. An optimal stopping problem with linear reward. Acta Mat., 132, 

111-151. 
[30] Wang, L., On the regularity theory of fully nonlinear parabolic equations L, Comm. 

Pure Appl. Math., 45 no. 1, 27-76, 1992. 
[31] Wang, L., On the regularity theory of fully nonlinear parabolic equations 11. , Comm. 

Pure Appl. Math., 45 no. 2, 141-178, 1992. 
[32] Wang, L., On the regularity theory of fully nonlinear parabolic equations HI., Comm. 

Pure Appl. Math., 45 no. 3, 255-262, 1992. 
[33] Wu, L.X., KwOK, Y. K., A front fixing finite difference method for the valuation of 

American options.. Journal of Financial Engineering, 6 (4), 83-97. 

Weierstrass Institute for Applied Analysis and Stochastics 

mohrenstr. 39 

10117, Berlin, Germany 

kampen@wias-berlin.de 



imsart-aap ver. 2007/12/10 file: aapamerican.tex date: July 7, 2008 



